Week 5 exercise: Simple MD simulations#
In this session we will implememnt simple MD simulations of increasing complexity.
Implement a numerical simulation of the dynamics of an Harmonic oscillator, using the verlet algorithm to integrate the equations of motion.
Using this notebook as a starting point (also accessible on Google colab):
Implement a dynamic simulation of an ideal gas in 2D.
Introduce an interaction \(\Gamma_{ii}(r)=r^{-12}\) potential to model a system of soft spheres, in a box defined by harminic repulsive walls.
Introduce a Lennard Jones interaction potential.
Observe what happens by performing simulations at different values of \(\epsilon\)
Introduce harmonic potentials to model diatomic molecules.
1D Harmonic Oscillator#
System Setup#
%matplotlib inline
## Import libraries to plot and do math
import matplotlib.pyplot as plt
import numpy as np
## Setup parameters of the HA
m=1.0; #mass
k=5; #Harmonic Constant
req=1.5; #Equilirbium
##
## Set the initial Conditions
r0=1.0; # Initial position
v0=.5; # Initial velocity
## Define the timestep and the total time
dt=0.01; #timestep
total_time=50;
# Compute the total number of steps
nsteps=int(total_time/dt); # Total number of steps
## Initialise vectors
v=np.zeros(nsteps)
r=np.zeros(nsteps)
time=np.zeros(nsteps)
POT=np.zeros(nsteps)
KIN=np.zeros(nsteps)
Definition of useful functions.#
Calculation of the force and potential energy#
The force is computed as the negative of the potential derivative with respect to the degree of the degree of freedom.
The potential is \(V(r(t))=\frac{1}{2}k(r(t)-req)^2\), hence the force is:
In the code below the function computing the force is implemented as:
force=lambda r, k, req, : -k*(r-req)
and the function computing the potential energy is defined as:
energy_pot=lambda r, k, req: 0.5*k*np.power(r-req,2)
Calculation of the velocity and kinetic energy#
The velocity at time \(t\), \(v(t)\) is computed as:
velocity=lambda r, r_past, dt: (r-r_past)/2/dt
The kinetic energy \(K(v)=\frac{1}{2}mv^2\)
energy_kin=lambda v, m: 0.5*m*np.power(v,2)
Propagation of the system’s dynamics:#
We use the Verlet algorithm:
implemented as:
verlet=lambda r, r_past, force, mass, dt: 2*r-r_past+(dt**2)*force/mass
## Useful functions
force=lambda r, k, req, : -k*(r-req)
energy_pot=lambda r, k, req: 0.5*k*np.power(r-req,2)
velocity=lambda r, r_past, dt: (r-r_past)/2/dt
energy_kin=lambda v, m: 0.5*m*np.power(v,2)
verlet=lambda r, r_past, force, mass, dt: 2*r-r_past+(dt**2)*force/mass
# initialise position and velocity
time[0]=0;
time[1]=dt;
v[0]=v0;
r[0]=r0;
r[1]=r[0]+v[0]*dt
v[1]=velocity(r[0],r[1],dt);
for i in np.arange(0,2):
POT[i]=energy_pot(r[i],k,req)
KIN[i]=energy_kin(v[i],m)
Compute Trajectory#
Propagate the dynamics with the Verlet algorithm :
# Compute trajectory
for ts in np.arange(1,nsteps-1): #Cycle over timesteps
f=force(r[ts],k,req) # Compute the force
r[ts+1]=verlet(r[ts],r[ts-1],f,m,dt) # Compute the next position
time[ts+1]=time[ts]+dt # update the clock
v[ts]=velocity(r[ts-1],r[ts+1],dt) # compute the velocity
POT[ts]=energy_pot(r[ts],k,req) # compute the Potential Energy
KIN[ts]=energy_kin(v[ts],m) # compute the Kinetic Energy
v[ts+1]=v[ts]
plt.plot(time,r,label='position');
plt.plot(time,v,label='velocity');
plt.legend();
plt.xlabel('time');
# Conservation of energy
plt.plot(time,POT+KIN,label='Total Energy')
plt.plot(time,KIN,label='Kinetic Energy')
plt.plot(time,POT,label='Potential Energy')
plt.ylim([0, np.max(POT+KIN)+0.05])
plt.legend()
plt.xlabel('time');
# Motion in phase space
plt.plot(r[10:],v[10:])
plt.xlabel('position');
plt.ylabel('velocity');
%%capture
# Visualize trajectory
%matplotlib inline
from matplotlib import animation, rc
from IPython.display import HTML
# Animate the results
def init():
line.set_data([], [])
return (line,)
fig, ax = plt.subplots(figsize=(8, 5))
#ax.set_xlim(( 0, total_time))
ax.set_xlim(-(np.max(r)/2+0.1*np.max(r)), np.max(r)/2+0.1*np.max(r))
line, = ax.plot([], [], lw=2, marker='o', markersize=45, markerfacecolor=(0.8, 1.0, 0.8, 0.5),
markeredgewidth=1, markeredgecolor=(0, 0, 0, .5), linestyle='-.',color='r')
def animate(i):
x = np.array([-r[i]/2,r[i]/2])
y = np.array([0,0])
line.set_data(x, y)
return (line,)
# call the animator. blit=True means only re-draw the parts that have changed.
anim = animation.FuncAnimation(fig, animate, init_func=init,
frames=np.arange(1,nsteps,50), interval=100, blit=True)
HTML(anim.to_html5_video())
---------------------------------------------------------------------------
RuntimeError Traceback (most recent call last)
Cell In[7], line 1
----> 1 HTML(anim.to_html5_video())
File ~/anaconda3/lib/python3.11/site-packages/matplotlib/animation.py:1285, in Animation.to_html5_video(self, embed_limit)
1282 path = Path(tmpdir, "temp.m4v")
1283 # We create a writer manually so that we can get the
1284 # appropriate size for the tag
-> 1285 Writer = writers[mpl.rcParams['animation.writer']]
1286 writer = Writer(codec='h264',
1287 bitrate=mpl.rcParams['animation.bitrate'],
1288 fps=1000. / self._interval)
1289 self.save(str(path), writer=writer)
File ~/anaconda3/lib/python3.11/site-packages/matplotlib/animation.py:148, in MovieWriterRegistry.__getitem__(self, name)
146 if self.is_available(name):
147 return self._registered[name]
--> 148 raise RuntimeError(f"Requested MovieWriter ({name}) not available")
RuntimeError: Requested MovieWriter (ffmpeg) not available
Model of a Fluid Phase in a 2D box: Template#
Work with the function forceij to introduce different interaction potentials:
Ideal Gas:#
\(\Gamma_{ij}=0\)
Soft sphere:#
\(\Gamma_{ij}=r_{ij}^{-12}\)
Soft sphere + harmonic:#
\(\Gamma_{ij}=k_{ij}(r_{ij}-r_{eq})^2+r_{ij}^{-12}\)
Lennard Jones:#
\(\Gamma_{ij}=4\epsilon\left[r_{ij}^{-12}-r_{ij}^{-6}\right]\)
Define Functions#
## Useful functions
verlet=lambda r, r_past, force, mass, dt: 2*r-r_past+(dt**2)*force/mass
forcebox=lambda x, boxx,boxk: np.greater(np.abs(x),boxx)*(-boxk)*x
#Define the system's box
boxx=10 #x dimension of the simulation' box
boxy=10 #y dimension of the simulation' box
boxk=1 #k constant for harmonic repulsive force
#Number of particles
N=10
#mass of the particles
m=np.ones(N)
######## "Force Field" Parameters #######
HS=1; # Repulsive soft potential
k=25.0; # Harmonic oscillator constant
req=1; # Harmonic oscillator equilibrium distance
KAPPA=k*np.zeros([N,N])
epsilon=0;
##Use this function to implement different potentials
def forceij(xi,xj,yi,yj,HS,KAPPA,req,epsilon):
r=np.sqrt((xi-xj)**2+(yi-yj)**2); #Distance
#Ideal Gas
dVdr=0
#Repulsive Wall
# dVdr=-12*HS/(np.power(r,13))
#Repulsive Wall + Harmonic potential
#... ... ...
#Lennard Jones Potential
#... .... .... ... ....
cx=-dVdr*((xi-xj))/r; #Pairwise component of the force in x
cy=-dVdr*((yi-yj))/r; #Pairwise component of the force in y
return [cx,cy]
def print_progress(iteration, total, bar_length=50):
progress = (iteration / total)
arrow = '*' * int(round(bar_length * progress))
spaces = ' ' * (bar_length - len(arrow))
print(f'\r|{arrow}{spaces}| {int(progress * 100)}% | ', end='', flush=True)
System Setup#
## Set the initial Conditions
# Random initial positions
x0=(np.random.rand(N)*2*boxx)-(boxx); #Initial position in x
y0=(np.random.rand(N)*2*boxy)-(boxy); #Initial position in y
# Random initial velocities
v0=(np.random.rand(2,N)-0.5); # Initial random velocitites
## Define the timestep and the total time
dt=0.005; # Timestep
total_time=100; # Total simulation time
nsteps=int(total_time/dt); # Total number of steps
## Initialise vectors
time=np.zeros(nsteps)
## Compute a trajectory with the Verlet Algorithm
# Initialise positions at t-dt
xp=x0;
yp=y0;
# Position at time t
x=xp+v0[0,:]*dt;
y=yp+v0[1,:]*dt;
# Position at time t+dt
xnew=np.zeros(N);
ynew=np.zeros(N);
# time
time=np.arange(0,nsteps);
time[0]=0;
time[1]=time[0]+dt;
## Initialize verctors for plotting
xx=np.zeros((np.size(time),N));xx[0]=x0
yy=np.zeros((np.size(time),N));yy[0]=y0
## |------------------|
## |Compute trajectory|
## |------------------|
for timestep in np.arange(1,nsteps): #Cycle over timesteps
timestep=int(timestep) #Make sure timestep is an integer
# Initialise force vectors
fx=np.zeros(N);
fy=np.zeros(N);
# Cycle over all particles
for i in np.arange(0,N):
fx[i]+=forcebox(x[i],boxx,boxk)
fy[i]+=forcebox(y[i],boxy,boxk)
for j in np.arange(i+1,N):
[cx,cy]=forceij(x[i],x[j],y[i],y[j],HS,KAPPA,req,epsilon)
fx[i]=fx[i]+cx; #update total x-component of the force on particle i
fx[j]=fx[j]-cx; #update total x-component of the force on particle j
fy[i]=fy[i]+cy; #update total y-component of the force on particle i
fy[j]=fy[j]-cy; #update total y-component of the force on particle j
xnew[i]=verlet(x[i],xp[i],fx[i],m[i],dt) # new position (x-component)
ynew[i]=verlet(y[i],yp[i],fy[i],m[i],dt); # new position (y-component)
print_progress(timestep,nsteps)
# Reassign positions
xp=x; yp=y; x=xnew+1-1; y=ynew+1-1;
## Store trajectory for animation
xx[timestep]=x;
yy[timestep]=y;
|**************************************************| 99% |
Visualization#
%%capture
## Display the trajectory
%matplotlib inline
from matplotlib.animation import FuncAnimation
from matplotlib import animation, rc
from IPython.display import HTML
fig, ax = plt.subplots(figsize=(10, 10))
line, = ax.plot([])
ax.set_xlim(-boxx, boxx)
ax.set_ylim(-boxy, boxy)
line, = ax.plot([], [], lw=2, marker='o', markersize=40, markerfacecolor=(0.8, 1.0, 0.8, 0.5),
markeredgewidth=1, markeredgecolor=(0, 0, 0, .5), linestyle=' ',color='red')
# initialization function: plot the background of each frame
def init():
line.set_data([], [])
return (line,)
def animate(frame_num):
x=xx[frame_num,:]
y=yy[frame_num,:]
line.set_data((x, y))
return (line,)
# call the animator. blit=True means only re-draw the parts that have changed.
anim = animation.FuncAnimation(fig, animate, init_func=init,
frames=np.arange(1,int(nsteps),50), interval=50);
HTML(anim.to_jshtml())